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Although the topic of fiducial marker-based alignment in electron tomography (ET) has been widely dis- 
cussed for decades, alignment without human intervention remains a difficult problem. Specifically, the 
emergence of subtomogram averaging has increased the demand for batch processing during tomo- 
graphic reconstruction; fully automatic fiducial marker-based alignment is the main technique in this 
process. However, the lack of an accurate method for detecting and tracking fiducial markers precludes 
fully automatic alignment. In this paper, we present a novel, fully automatic alignment scheme for ET. 


es jah Our scheme has two main contributions: First, we present a series of algorithms to ensure a high recog- 
Alignment erapny nition rate and precise localization during the detection of fiducial markers. Our proposed solution 


reduces fiducial marker detection to a sampling and classification problem and further introduces an 
algorithm to solve the parameter dependence of marker diameter and marker number. Second, we pro- 
pose a novel algorithm to solve the tracking of fiducial markers by reducing the tracking problem to an 
incomplete point set registration problem. Because a global optimization of a point set registration 
occurs, the result of our tracking is independent of the initial image position in the tilt series, allowing 
for the robust tracking of fiducial markers without pre-alignment. The experimental results indicate that 
our method can achieve an accurate tracking, almost identical to the current best one in IMOD with half 
automatic scheme. Furthermore, our scheme is fully automatic, depends on fewer parameters (only 
requires a gross value of the marker diameter) and does not require any manual interaction, providing 
the possibility of automatic batch processing of electron tomographic reconstruction. 

© 2015 Elsevier Inc. All rights reserved. 


Point set registration 
Four point invariant 


1. Introduction collection. Therefore, to obtain high-quality three dimensional 


density maps, the projection parameters of the tilt series should 


Electron tomography (ET) techniques have become an indis- 
pensable tool in structural biology. In ET, the three-dimensional 
density of the ultrastructure is reconstructed from a series of 
micrographs (tilt series) taken in different orientations. Generally, 
the geometrical information about the tilt series can be recorded 
by the electron microscope. However, the projection environment 
may be affected by mechanical instability, and inevitable transfor- 
mations and deformations of the sample occur during data 
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be calibrated accurately before reconstruction. 

Two main types of alignment methods are available for ET: 
marker-free alignment and marker-based alignment. Each type of 
method has its own applicability and limitations. Marker-free 
alignment (e.g., cross-correlation (Guckenberger, 1982) and com- 
mon lines (Liu et al., 1995) used in coarse alignment, the iterative 
alignment methods combining cross-correlation with reconstruc- 
tion and reprojection (Winkler and Taylor, 2006, 2013), and 
feature-based alignment methods that replace fiducial markers 
with features (Brandt et al., 2001a; Brandt and Ziese, 2006; 
Castafio-Diez et al., 2007, 2010; Phan et al., 2009; Sorzano et al., 
2009; Han et al., 2014)) does not require fiducial markers 
to be embedded in the sample. However, the applications of 
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marker-free alignment are still limited by Signal-to-Noise Ratio 
(SNR), the intrinsic biological structure of the sample and the con- 
trast of the micrographs. Fiducial marker-based alignment 
(Lawrence, 1992; Kremer et al., 1996; Frank, 2006) requires that 
fiducial markers be embedded in the sample. Because fiducial 
markers have a high contrast with the background, the 
positions of the markers can be determined precisely, which 
further improves the alignment accuracy. Therefore, fiducial 
marker-based alignment is most accurate and has been widely 
used for high resolution electron tomography. Particularly for 
cryo-ET datasets with low SNR, fiducial marker-based 
alignment could be the only choice (Amat et al., 2008). However, 
to be noted that, the embedding of fiducial markers may interfere 
with the sample and introduce undesirable artifacts in the 
reconstruction. 

Fiducial marker-based alignment requires the positions and the 
corresponding relationship of markers in different tilted micro- 
graphs. Thus, marker-based alignment should first localize the 
embedded fiducial markers in micrographs and then track these 
markers throughout the tilt series. Although fiducial marker- 
based alignment has existed almost as long as structural biology, 
there are still no robust automatic schemes to perform the 
detection and tracking of fiducial markers. With the advent of 
the subtomogram averaging technique (Briggs, 2013), batch 
tomographic reconstruction becomes an urgent demand, and fully 
automatic alignment is the main technique for batch tomographic 
reconstruction. 

The manual intervention in marker detection and tracking is 
the bottleneck of marker-based alignment for high-throughput 
tomographic reconstructions. A number of studies have focused 
on the marker detection and tracking issue. The fiducial mar- 
ker detection is usually performed by template matching 
(Kremer et al, 1996; Brandt et al., 2001b; Zheng et al., 
2007; Amat et al., 2008; Cao et al., 2011). However, the result 
of template matching is usually affected by the value of the 
threshold used to determine the similarity between the image 
patches and the fiducial marker template. Furthermore, the 
value of the threshold used in the template matching usually 
varies for different datasets, and an unsuitable threshold 
results in an unacceptable output. Current studies on fiducial 
marker-based alignment also have invested minimal effort in 
the refinement of the marker diameter parameter. Usually, 
the diameter parameter is set by users, and an incorrect value 
also can corrupt the result. The fiducial marker tracking is usu- 
ally performed by local search between consecutive images, 
using the underlying local geometry information or local 
similarity of patches (Ress et al., 1999; Brandt et al., 2001b; 
Amat et al., 2008; Cao et al., 2011). The reason why previous 
studies employ local search to track fiducial markers is that if 
the problem is modeled as a globally consistent tracking, the 
required computational resources will be extremely large. The 
use of local geometry information can improve tracking speed, 
but it cannot guarantee global consistency. Here, we present a 
novel fully automatic scheme for fiducial marker-based 
alignment. 

First, we developed a complete algorithm to solve the detec- 
tion problem automatically. Two types of parameters are used 
in the detection of fiducial markers. The first type includes the 
information of fiducial markers, such as the marker diameter 
and the shape of the fiducial marker. The second type of param- 
eters includes the thresholds used to detect of fiducial markers. 
A robust automatic detection of fiducial markers depends on 
reasonable adjustment and estimation of these parameters. First, 
we developed an algorithm to refine the marker diameter. 
Then, we convert the maker detection problem into a sampling 
and classification problem. By solving this problem, we can 


automatically detect the fiducial markers, estimate the number 
of markers in each micrograph and eliminate the dependence 
of the threshold used in the detection. Thereafter, we developed 
an algorithm to refine the localization precision of the detected 
markers. 

After fiducial marker detection, we propose a global 
optimization algorithm to solve the tracking problem. Under 
the weak perspective projection model, corresponding image 
features in two respective 2D views from the identical planar 
surface are related by an affine transformation (Huttenlocher 
and Ullman, 1990; Koenderink and van Doorn, 1991). Based 
on the affine transformation relationship of the two projections, 
we reduce the tracking problem to an incomplete point set reg- 
istration problem and propose a novel method to solve the 
matching of corresponding fiducial markers. Unlike traditional 
studies (Ress et al., 1999; Brandt et al., 2001b; Amat et al., 
2008; Cao et al., 2011) that treat the matching of corresponding 
fiducial markers as a tracking problem and use the local geom- 
etry or similarity of patches to accelerate the tracking, our solu- 
tion does not use the local geometry information. Because the 
optimization of the point set registration is globally consistent, 
our tracking method is independent of the initial state of the tilt 
series and allows inputs with unrefined fiducial marker 
positions (possibly contaminated by outliers and missing data). 
Furthermore, the optimization of the point set registration can 
be solved by common computational resources in an acceptable 
time. 

Finally, we combine the detection and tracking methods 
into a fully automatic alignment scheme and evaluate the 
scheme with experimental data. The theoretical analysis and 
experimental results show the efficiency of the proposed 
scheme. 


2. Method 


The steps of our alignment scheme are illustrated in Fig. 1.! 
The first step is to refine the parameters of the fiducial markers 
and detect the fiducial markers in the tilt series. We develop an 
algorithm to estimate the diameter parameter of the fiducial 
markers. A novel sampling and classification algorithm is used 
to ensure the exhaustive detection of the fiducial markers, and 
a new algorithm used to refine the final positions of the extracted 
fiducial markers is also illustrated as a post-detection subprocess. 
The second step is to match corresponding markers. We design a 
novel algorithm based on the idea of incomplete point set regis- 
tration to guarantee a globally consistent matching of two marker 
sets. The third step is to track matching pairs consistently across 
the tilt series and to stretch the tracks as long as possible. A 
strategy that accelerates the tracking and compensates for the 
missing of fiducial markers is used in this step. The final step 
is to optimize the projection parameters with the determined 
tracks and, when necessary, to transform the images 
geometrically. 


2.1. Marker parameters refinement and marker detection 


To detect fiducial markers, we need a clear definition of the 
appearance of fiducial markers in ET: the shape of a fiducial mar- 
ker in a micrograph is a circular or nearly circular area with an 
outstanding contrast to the nearby area and the shape of an 
identical fiducial marker does not change substantially with 


1 Our system uses the inverted color image. For example, if the pixel value ranges 
from 0-1, 0 presents white and 1 presents black; if the pixel value ranges from 0-255, 
0 presents white and 255 presents black. 
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Refine fiducial marker parameters 
& detect fiducial markers 


Match sets of 
fiducial markers 


Track corresponding markers 
& extend tracks 


Optimize projection parameters 


Align stack 


Fig. 1. Flowchart of our automatic alignment scheme. 


different tilt angles. In addition, we notice that the theoretical 
capacity, i.e., the number of fiducial markers that a micrograph 
is able to contain, is limited by the image size and the diameter 
of the markers (regardless of overlapping). In traditional studies 
(Ress et al., 1999; Brandt et al., 2001b; Amat et al., 2008; Cao 
et al., 2011), the number of fiducial markers is set by the users, 
and the candidates of fiducial markers are ranked according to 
their Normalized Cross-Correlation (NCC) value. When the rank 
of a candidate is lower than the number inputted by the user, 
the candidate is assumed to be a fiducial marker. This type of 
strategy works well in most conditions but is not robust against 
incorrect user inputs. Furthermore, when the user has input a 
large number (for example, 10,000 fiducial markers), no 
applicable threshold develops. A hard threshold used to deter- 
mine the final set of fiducial markers also is not applicable 
because of the change of amplitude of micrographs in different 
tilt angles. 

Before further discussing our method, we would like to 
emphasize the theoretical capacity of fiducial markers in a 
micrograph. Assuming we already know the exact value of the 
marker diameter (denoted as d), given a micrograph with width 
w and height h, the theoretical maximum number of fiducial 
markers that can be detected from the micrograph is 


(w-h)/@. Furthermore, if we have already picked out all the 


(w-h)/d? candidate patches from the micrograph, our work is 
to classify these patches and determine the real fiducial mark- 
ers. Therefore, in the case that the exact value of the marker 
diameter is already known, we can reduce the detection prob- 
lem to a classification problem without considering the actual 
number of fiducial markers or the threshold. On the contrary, 
if we have a method to classify the candidate patches and have 
already extracted all the fiducial markers, we can determine the 


most accurate value of marker diameter. Therefore, we propose 
a novel algorithm to detect the fiducial markers. The detection 
algorithm consists of two sub-algorithms: a fiducial marker 
detecting algorithm based on the presupposition of known mar- 
ker diameter and a marker diameter refining algorithm based 
on the presupposition of known exact fiducial markers. We will 
first introduce the marker detecting algorithm (Algorithm 1) 
and then introduce the marker diameter refining algorithm 
(Algorithm 2). 


Algorithm 1. Detect fiducial markers under a given value of 
marker diameter 


procedure: DetectFipuctaLManxers(d,img) 

: Normalize the amplitude of img 

: Create fiducial marker template tmplt 

: corr — NCC(img, tmplt) 

: Normalize the amplitude of corr 

: CANDI — Ø 

: thre1 — AVG(corr), thre2 — AVG(img) 

: for i = 0 to w in step length of d do 

for j = 0 to h in step length of d do 

9: Find the point p;;) with maximum NCC value 

10: if corr(i,j) > thre1 and img(i,j) > thre2 then 

11: Paj- ncc — corr(i,j), paj val — img(i,j) 

12: CANDI — pj) U CANDI 

13: end if 

14: end for 

15: end for 

16: Analysis the distribution of ncc and val of candidates {In 
our implement, we fit a Gaussian distribution model} 

17: Calculate the mean u and the standard deviation o of 
{Piij)-NCC x Paj val} 

18: thre3 — u + 2.50 

19: FID — Ø 

20: for all p;; € CANDI do 

21: if paj-ncc x paj-val > thre3 then 

22: FID — FID U Pij) 

23: end if 

24: end for 

25: Test the shape of fiducial markers in FID by region 
growing 

26: Refine the position of fiducial markers in FID 

27: return FID 


ANnAUBRWN 


The flowchart of the marker detecting algorithm is illustrated 
in Fig. A.1. The marker detecting algorithm requires two inputs: 
the micrograph and a template of fiducial markers. The template 
can be an average of fiducial markers or an artificial template 
generated from the already known marker diameter. As discussed 
above, the fiducial markers have very different shapes and 
contrasts strongly with the background. Here, the circular shape 
and the contrast is used as two characteristics of the fiducial 
markers in our solution. The algorithm for a given diameter (arti- 
ficial template is used as marker template) is described in detail 
in Algorithm 1, where d is the input value of the diameter (unit: 
pixel), img is the original micrograph, FID is the set of fiducial 
marker positions, AVG(-) is the operation to obtain the mean pixel 
value of an image and the set CANDI represents the candidate 
points of fiducial markers. First, an artificial template is created 
from the input diameter information (a black circular area with 
a white background). Then, the Normalized Cross-Correlation is 
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Fig. 2. The distribution of the shape labels and contrast labels of candidate fiducial markers. 


used to determine the similarity between patches of the micro- 
graph and the template. The amplitudes of img and corr are nor- 
malized to 0 ~ 1 for further analysis (without cutoff).? The image 


is then divided into (w-h)/d’ patches and for every patch, the 
point with the maximum cross-correlation value is selected as a 
candidate of fiducial marker position. The selected marker positions 
whose pixel value in micrograph is smaller than the mean pixel 
value of img or corresponding NCC value with the template is smal- 
ler than the mean value of corr will be excluded from the candidate 
set. In other words, a very loose threshold is used to further refine 
the candidate set. The NCC value of the patch is assigned to the 
candidate as the ‘shape label’ (denoted as pj;).ncc) and the 
corresponding pixel value of the point in the micrograph is 
assigned to the candidate as ‘contrast label’ (denoted as p; val). 
It should be noted that in set CANDI, there are still thousands can- 
didate points remaining after being filtered by the loose threshold. 
We use the candidates extracted from the zero-tilt slice of the Cen- 
triole (d = 8 pixels), an example downloaded from the website of 
IMOD (Kremer et al., 1996), as our illustration to show the charac- 
teristics of fiducial markers. Fig. 2 illustrates the typical distribution 
of the values of the shape label and contrast label of these candi- 
dates. As shown in Fig. 2(a), the background and the real fiducial 
markers can be easily separated by a boundary line. However, to 
make the computer understand the difference between the back- 
ground and real fiducial markers, we still need other technique 
supports. We have tried cluster analysis, for example, k-means 
and OPTICS (Ankerst et al., 1999), to separate the background and 
real fiducial markers. Though these techniques produce good 
separation and figure out the final candidate of fiducial markers, 
they cost relatively high computational resources. Our used solu- 
tion is a threshold-based method relying on the statistical proper- 
ties of the ‘labels’, which can reach comparable accuracy in a 
simple way. In this solution, for every candidate, a score 
s(p) = NCC valuex pixel value is calculated. Given a threshold 
whose value is u+ 2.50, where u is the mean of s(p) and a is 
the standard deviation of s(p), if s(p) is larger than the threshold, 
the corresponding candidate will be reserved. Fig. 2(b) illustrates 
the histogram of the NCC value x pixel value for all candidates. It 
can be found that the distribution of the values of s(p) approxima- 
tively follows a Gaussian distribution, which is in agreement with 
the central limit theorem. Maybe some good candidates are 
excluded by this solution, however, the proportion is very small 


2 For the input image, the details of normalization is carried out as following steps: 
1. Calculate the average value (denote as avg) and the standard deviation (denote as 
std) of the pixels of the image. 2. Get the maximum pixel value (denote as max) and 
minimum pixel value (denote as min) of the image. 3. Calculate the bound value for 
normalization: low bound low=MAX(avg—3-std,min); upper bound 
upper = MIN(avg + 3 - std, max). 4. For each pixel I(x,y), the value is recalculated as 
I(x, y) = (I(x, y) — low) /(upper — low). 


and they can be recovered in marker tracking and extension. 
Finally, region growing is used as the final check to exclude false 
and low quality markers? and the positions of fiducial markers 
are recalculated (the refinement of fiducial marker positions con- 
sists of two steps: the first step is to relocalize the fiducial markers 
by centroid calculation; the second step is to refine the positions by 
our position refinement algorithm, i.e. Algorithm 3, which will be 
discussed below). 


Algorithm 2. Refine the diameter of fiducial markers 


procedure: ReFINEMArKERDIAMETER(d,img) 

1: while ¿ > éo do 

2: Initialize search space of diameter 
DIA = {di|d — €/2 < dj < d+ €/2,i=1,2...M} 

3: for all d; € DIA do 

A: Create fiducial marker template tmplt 

5: FID = DetectFiducialMarkers(d;, img) 

6: Create the average of fiducial markers avg_tmplt from 
FID with diameter d;i 

7: d;.sim — NCC(avg_tmplt, tmplt) 

8: end for 

9: Find the diameter with maximum score dmax and sub- 
max score d ub 

10: de (dmax + dsup)/2, č — (dmax — dsub) 

11: end while 

12: return d 


Although the diameter of fiducial markers can be determined 
using edge detection methods (for example, the Sobel Filter 
Kremer et al. 1996), the performance of edge detection is 
negatively affected by the SNR of cryo-ET datasets. We propose 
an algorithm to guarantee a robust estimation of the marker 
diameter. The flowchart of one loop of the algorithm is illustrated 
in Fig. A.2, and the detailed algorithm is described in Algorithm 2, 


3 Region growing is a simple pixel-based image segmentation method, which can 
identify image region with close values of pixels. Region growing needs two 
parameters: the initial seed points and the threshold to determine whether the 
values of two pixels is close. Here, the region growing is used to get the region of 
fiducial markers from the original micrograph. For each candidate, the center position 
(denoted by p,;;)) of marker candidate is inputted as the initial seed point. Assuming 
the given candidate of fiducial marker is a good one, the whole fiducial marker will be 
limited in the image patch (d x d size) centered in p; j. Because we have extracted the 
center positions of the candidate fiducial markers and already known the value of 
marker diameter, for each candidate, the threshold used in region growing can be 
determined by the statistics of pixel values in its corresponding image patch (d x d 
size). For a candidate fiducial marker, if the width and height of the region finally got 
by region growing is too small or too large compared with the already known 
diameter value, the candidate will not be considered. Therefore, the region growing in 
fact rechecks the contrast of candidate fiducial markers to the background and the 
diameter of candidate fiducial markers. 
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where d is the input value of the diameter, img is the original 
micrograph, é is the search range and é is the exit condition. 
In our implementation, a grid search is used in the algorithm, 
where d; is generated by an internal step that will change in each 
iteration and M is the size of discretized search space DIA. First, 
the value of the diameter obtained from the search space is 
assumed to be correct, and the corresponding fiducial markers 
are detected by Algorithm 1. Then, all of the fiducial markers in 
set FID are assumed to be real fiducial markers and localized 
accurately; therefore, an average template of fiducial markers 
can be obtained by averaging all the fiducial markers. A similarity 
score is obtained by calculating the NCC between the artificial 
template tmplt and the average template avg_tmplt. The above 
process is repeated until all the similarity scores for the diameters 
in the search space are calculated. The search space is updated by 
the diameters with the maximum value of similarity score. 
Finally, when the search range is smaller than ¢éo, the iteration 
stops and an estimation of the diameter is outputted. Typically, 
the zero-tilt slice of the micrographs is used as a reference for 
the diameter refinement. Our experiments indicate that an čo 
value of 0.5 pixels is sufficient. 

By introducing Algorithms 1 and 2, we have solved the problem 
of the dependence of parameters of marker diameter and the 
number of fiducial markers. When a suitable diameter value is 
determined in Algorithm 2, a fiducial marker template (average 
template) can be generated by averaging the detected fiducial 
marker patches; this average template can replace the artificial 
template used in Algorithm 1 and further improve the performance 
of our method. An illustration of the performance of our algorithm 
will be provided in the Section 3. 


Algorithm 3. Refine the positions of fiducial markers 


procedure: RerineEMarkerLocaization(d,img, FID) 

1: res — int_max {residual} 

2: while res > resin do 

3: Get all the patches that represent fiducial markers 
(image area with d x d size, centering in the position of 
markers) from FID, denote the patches as FID IMG 

4: Enlarge the patches in FIDIMG K times (K = 16, in our 
implementation) 

5: Create the average of fiducial markers avg_tmplt from 
FID_IMG (Because FID_IMG have been enlarged, avg_tmplt is 
also enlarged K times) 


6: res~—O 

7: forall fid; € FID do 

8: Get an image area of a marker ((d + 2 - B) x (d + 2 - B) 
size, centering in the position of fid;), denote the image area 
as patch; 

9: Enlarge the patch; K times 

10: corr — NCC(patch;,, avg_tmplt) 

11: Find the point p; with maximum NCC value in corr 

12: map p; to the position in img, update the value of p; 

13: res — res + DISTANCE(fid,, pi) 

14: fid; — p; 

15: end for 


16: res — res/(size of FID) 
17: end while 
18: return FID 


Localization precision is another issue. Although the centroid of 
fiducial markers can provide sub-pixel marker localization, its 
performance is easily affected by the background of micrographs 
(the centroid only considers the pixel value). Here, we propose 
an algorithm to refine the localization of fiducial markers by 


iteratively re-localizing the fiducial markers, generating an average 
template and re-scoring the similarity. The flowchart of the 
algorithm is illustrated in Fig. A.3 and the detailed algorithm is 
described in Algorithm 3. In our opinion, the shape is more impor- 
tant than the contrast for the final localization of fiducial markers. 
First, we generate an average template of fiducial markers that 
have been enlarged K times (1/K represents the measurement 
precision; K=16 in our implementation). The image areas 
(patch; corresponds to fid;) around a marker position are also 
selected and enlarged. We will spare some boundary (size in B) 
for further adjustment of the marker localization. Based on the 
NCC value of each patch, the position of each fiducial marker fid; 
is refreshed (p,;, in minimum scale 1/K) and the average residual 
res between all the fid; and p; is calculated. This procedure is per- 
formed iteratively until the value of res is smaller than resin 
(typically we set reSmin to 1/K). Our algorithm can localize the 
position of fiducial markers to a minimum scale of 1/K, and our 
algorithm sufficiently considers the shape of the fiducial markers 
during the re-localization. Furthermore, because the average tem- 
plate of fiducial markers is recalculated for each micrograph, our 
algorithm can suppress the variations in the scale and contrast of 
fiducial markers caused by the change of tilt angle. 


2.2. Corresponding fiducial marker matching 


Two characteristics should be considered in fiducial marker 
matching. The first consideration is that the positions of fiducial 
markers in micrographs can change substantially because of the 
perspective projection and in-plane rotation of samples. The sec- 
ond consideration is that there may be missing or falsely detected 
fiducial markers for the micrographs from different tilt angles. 

First, we restate the problem as follows: Let ‘point set’ denote 
the positions of the fiducial markers extracted from a projection. 
Given two point sets (M and S) belonging to different views, an 
affine transformation T(-) applied to the moving “model” point 
set M can be found so that there exists a subset of T(M) with 
the maximum cardinality in which the points are corresponding 
to the points from a subset of the static “scene” set S under a select 
measure of distance. The corresponding relationship of these 
points is called ‘congruent’ and the two subsets with maximum 
corresponding points to each other are called maximum congruent 
subsets. The subset of T(M) with maximum corresponding peers is 
denoted as largest common point set (LCP). 

Because the cardinalities of M and S are not equal (M and S 
represent the detected fiducial marker positions of different micro- 
graphs and the size of detected fiducial markers from different 
micrographs are not equal), it is an incomplete point set registra- 
tion problem for fiducial marker tracking. Although there are many 
works discussing the problem of point set registration (Bishnu 
et al., 2006; Rusinkiewicz and Levoy, 2001), most of the works 
focus on the case of rigid transformation. Recently, some efforts 
based on heuristic methods (Zhang et al., 2003) have been con- 
ducted to solve the affine point set registration problem. However, 
the performance of these methods heavily depends on the choice 
of parameters and is limited by the application field. To solve the 
incomplete point set registration problem, we must make a trade- 
off between the robustness and computational complexity. 

Before further analysis of our solution, we introduce three basic 
concepts. The first concept is the affine invariant ratio of four 
points (Huttenlocher, 1991), which has been widely used in com- 
puter graphics and vision (Aiger et al., 2008). As shown in Fig. 3, 
let points a, b, c, d compose a point set S4. Point e is the intersection 
of line ab and cd. If there is a point set S2 generated from Sı under 
the transformation T(-), then the ratio rı = ||a —e||/||a — b|| and 
t2 = ||c —e||/||c—d|| are preserved under affine transformation, 
i.e., |æ —e'||/\la’ —b'|| = rı and |e — e'||/||c’ — d'|| = r2. The use of 
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d? 


Fig. 3. Illustration of a four-point congruent set. 


the affine invariant can significantly reduce the computational 
complexity. The second concept is the random sample consensus 
(RANSAC) (Fischler and Bolles, 1981). RANSAC is a widely used 
technique for the robust estimation of data with noise and outliers. 
Our matching algorithm is based on the idea of RANSAC. The third 
concept is the affine transformation model. Micrographs recorded 
in a tilt series obey the laws of affine projection (affine projection 
is a special case of weak perspective projection), the relationship of 
corresponding image points can be modeled as an affine transfor- 
mation (Huttenlocher and Ullman, 1990; Koenderink and van 
Doorn, 1991). Transformation T(-) is modeled as x’ = Ax + t, where 
x and x’ are vectors representing points in R?, A is a 2 x 2 matrix, 
and t is the translation. 


Algorithm 4. Find the maximum congruent subsets of M and 
S under appropriate measure of distance 


procedure: FinpMatcu(M, S, dist_thre) 

1: Generate all the quadrilaterals Q; with wide-base, 
® = {Q,|Q; is a 4 point subset ofS} 

2: Calculate invariant ratios of all Q;, i.e. Q;.inval and 
Q;.inva2 

3: Create kd-tree of ® with Q;.inval and Q;.inva2 as search 
keys, denote the kd-tree as KD(®) 

4: Calculate the maximum iteration count max_count 

5: iteration — 0 

6: while iteration < max_count do 

7: Randomly select a quadrilateral Pyan from M 

8: Search KD(®), get all the Q; whose keys appropriately 
equal to Pran, denote the quadrilateral set as Y 

9: for all Q; € ¥ do 


10: Estimate the transformation T(-) to be applied to 
Pran, SO that T(P;an) appropriately equals Q; 

11: Apply T(-) to M 

12: Find C(T(M)) c T(M) and C(S) c S, so that C(T(M)) 


is congruent to C(S) with the maximum cardinality under 
distance threshold 3dist_thre 


13: if the cardinality of C(T(M)) is large enough then 

14: Re-estimate T(-) by C(M) and C(S), so that 
T(C(M)) appropriately equals C(S) 

15: Apply T(-) to M 

16: Find C(T(M)) c T(M) and C(S) c S, so that 


C(T(M)) is congruent to C(S) with the maximum cardinality 
under distance threshold dist_thre 


17: end if 

18: if C(T(M)) has the current maximum cardinality 
then 

19: MATCH_-PAIR — Pair < C(M),C(S) > 

20: Re-calculate the maximum iteration count 
max_count 

21: end if 

22: end for 

23: iteration — iteration + 1 


24: end while 
25: return MATCH _PAIR and T(-) 


Our solution is described in Algorithm 4, where M is the 
moving point set, S is the fixed point set, and dist_thre is the dis- 
tance threshold used to find the largest common point set. In 
our implementation, only the quadrilaterals with wide-base are 
used.* Although there are many quadrilaterals in a point set, most 
quadrilaterals are not suitable for estimating T(-). Matching with 
wide-base is more stable. Kd-tree is used to accelerate the search 
of quadrilaterals with appropriate keys (Bentley, 1975). The com- 
plexity of building a kd-tree differs according to its implementa- 
tion and generally refers to O(anlogn). The search complexity is 
O(logn) on average. Here, « is 2, referring to the two ratios of 
the 4-point affine invariant. The matching procedure based on 
RANSAC in our algorithm is as follows: First, a quadrilateral is 
random selected from M, a search is performed for the potential 
congruent point set of the selected quadrilateral in S, and the cor- 
responding transformation T(-) is calculated. Then the calculated 
transformation T(-) is applied to M, and the number of points 
in T(M) that is congruent with S is counted. If the cardinality is 
sufficiently large, accept the transformation T(-) as a solution. 
Repeat this process until the iteration reaches the maximum 
count. Typically, the transformation T(-) is estimated by least 
squares estimation. However, least squares estimation is not 
robust to noise. To overcome this disadvantage, a two-stage esti- 
mation of the transformation is used. In the first stage, T(-) is esti- 
mated using a randomly selected quadrilateral Pra, and its 
congruent set Q;. Then, a loose distance threshold 3 - dist_thre is 
used to count the appropriate congruent points in M and S. If 
the cardinality of the congruent subset C(T(M)) is sufficiently 
large, go to the second stage. In the second stage, the transforma- 
tion T(-) is re-estimated by C(M) and C(S), and the congruent sub- 
set is re-counted under a tight distance threshold dist_thre. In our 
solution, the maximum iteration count is refreshed according to 
the formula L = log(1 — p,)/log(1 — pł), where p, is the required 
success probability, p, is the probability that the randomly 
selected point from M is also present in S, and N is the 
cardinality of the randomly selected data set (N = 4). Let m be 
the number of points in M and n be the number of points in S. 
The generation of all quadrilaterals with wide-base is O(n? +k), 
where k is the total number of quadrilaterals with wide-base. 
Here the value of k depends on the diameter used as the 
wide-base. Typically, k is much smaller than n* if an appropriate 
value is selected for the diameter. Thus, the complexity of 
building kd-tree is O(klogk), which is smaller than O(n? logn). 
Considering that the complexity of selecting random quadrilaterals 
in M is O(m), the total complexity of our method is bound in 
max(O(n? + k),O(klogk),O(Lmlogn)). Because the p, in L is 
very large in our problem (the appearance of fiducial markers in 
different micrographs is stable), the value of L is relatively 
small. Therefore, our solution is much faster than an 
exhaustive method (approximately O(n) according to Zhang 
et al. (2003)). 


2.3. Fiducial marker tracking 


If only the corresponding relationship of markers in two views 
is found, the tracking of fiducial markers in different micrographs 
is easy. First, we employ a multi-level matching strategy presented 


4 A quadrilateral with wide-base means that the 4 points compose a convex 
quadrilateral and the diagonals of the quadrilateral are wide-base. Wide-base is a 
measure that judges from the distance of two points or the length of a line: Let us 
denote the maximum of the distances between two arbitrary points inside a shape by 
Dmax, if the distance of two points compared with Dmax is not very small, for example, 
> 0.4 - Dmax, it is said that the two points are wide-base. 
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in Han et al. (2014) to obtain a sufficient number of pairs of 
corresponding fiducial markers: 


(1) Initialize step = 1, while step <=MAXSTEP; 

(2) For the nth micrograph of the series, match the markers 
between the nth and (n + step)th point sets(n = 1,2...); 

(3) step = step + 1; 

(4) Repeat until the condition in (1) is not satisfied. 


This type of strategy has been widely used in ET (Brandt and Ziese, 
2006; Amat et al., 2008; Han et al., 2014). Experiments have shown 
that a value of 2 is acceptable for MAXSTEP in the situation of fidu- 
cial marker-based alignment.° After obtaining a sufficient number 
of pairs, we would like to use the tracking model (a data structure 
mainly copes with random insertion, searching of markers and 
matching collision) introduced in the work of Han et al. (2014) to 
accelerate the generation of tracks. The case in which a local pair 
of corresponding markers is missing is also easy to handle. Because 
we have already obtained the transformation T(-) of micrograph 
pairs during the marker matching, given an already detected fiducial 
marker, we can reasonably postulate the possible location of the 
missing marker in the next micrograph. Using template matching 
in the restrained area around the postulated location, we can easily 
find the missing fiducial markers. All potential missing markers are 
checked and re-detected in the tilt series until the tracks cannot be 
extended. It should be noted that for some markers in micrographs 
with high tilt angle, the tracks are naturally shorter than others 
because these markers have not been projected in the micrographs 
with low tilt angles. 


2.4. Projection parameter optimization 


Once all of the tracks have been generated, a suitable projection 
model is needed to determine the projection parameters. 
Generally, the photographing of ET micrograph is assumed to obey 
the law of affine projection. The projection model is formulated as 
follows (for convenience without the loss of generality, the 
coordinate and vector used here are given in homogeneous 
coordinates): 


m X 
R Y 
=K(1.0)( 5 D 7 (1) 
Ze 
1 


5 For MAXSTEP = 2, it is possible that tracking of some markers is inconsistent. For 
example, a marker reached by tracking from marker i on micrograph n to micrograph 
n+1and from n+ 1 to n + 2, need not be the same as the marker reached by tracking 
from marker i on micrograph n directly to micrograph n + 2 (step = 2). However, it is 
not a big problem, we have a solution that can solve ordinary cases. In our scheme, 
the matching of corresponding fiducial markers is carried out by point set 
registration. In point set registration, we also get the transformation T(-). When 
collision happens, we can use the obtained transformation T(-) to determine which 
one is better. Since the position of marker i, the transformation function from 
micrograph n to micrograph n + 1 (denote as Tn n+1) and the transformation function 
from micrograph n+1 to micrograph n+2 (denote as Ty1n.2) are known, the 
transformation function from micrograph n to micrograph n+2 (denote Tyn+1.n+2) 
can be calculated by Tnn+1 and Tn+1.n+2. Also, we have the transformation function 
directly from micrograph n to micrograph n+2 (denote T,,,,2). We can first use 
Trnsinz2 and marker i to get a predicted marker position p; and calculate the 
deviation from the measured marker position to p;; then we use Tp n+2 and marker i to 
generate a predicted marker position př and calculate the deviation from the 
measured marker position to p/. The position with smaller deviation to the predicted 
one is selected as the ‘real’ one. In some special cases, especially the measured two 
fiducial markers are very close to each other, our solution may not identify which is 
the real one. However, because most projection parameter estimation methods are 
robust to outliers, we tolerate a small number of fault case. 


where K is the intrinsic matrix of the camera, (I,0) is the projection 
matrix, R is an 3 x 3 orthogonal matrix, which can also be written as 
R,,- Rg - Ry, representing rotation in space, 


1 0 0 
R,=|0O cosa sing |, 
0 -sing cose 
cosp 0 -sinf 
Rp = Oo 1 0 ; 
sinf 0 cosß 
cosy siny 0 
R, = | -siny cosy 0 
0 0 1 


and t is the 3-dimensional translation vector of the sample. 
(X,Y,Z,1) is the spatial locations of the markers in the 
sample, and (u,v,z,) are the projections in the micrographs, i.e., 
(U/Ze, V/Z-) represents the positions of fiducial markers (x,y) in 
the micrograph. Generally, the intrinsic matrix K is fixed because 
the strength of electron beam and the intrinsic parameters of cam- 
era changes minimally during experiment. The following cost func- 
tion is used to optimize the parameters: 


B= (Proj(X;) = xij) ii 2) 
ij 


where {X j} is the set of spatial points to be estimated (in inhomoge- 
neous coordinates). Proj(-) is the process of projection defined in Eq. 
(1). x; is the configured position of fiducial marker and 6); 
(ôij € {0, 1}) indicates whether the jth marker is visible in micrograph 
i. Our aim is to find the most appropriate R and t with the objective of 
minimizing the cost function E. Generally, the initial values of the 
optimized parameters are assigned as follows: The rotation parame- 
ters «and y are set to zero, and £ is equal to the assumed tilt value. Like- 
wise, the position of the camera in the lowest absolute value of the tilt 
angle is chosen as the original coordinates. t is initialized according to 
the scene depth and initial tilt angle. K is determined based on the 
intended application. For each spatial point (X, Y, Z, 1), the initial value 
is calculated by triangulating projections in low-tilt views. 

Because the projection transformation and the spatial location 
of the markers are both unknown a priori and the transformation 
of the sample further affects the projection location, the optimiza- 
tion of E is a non-linear problem. Many mathematical techniques 
are available to solve the non-linear problem (Lawrence et al., 
2006; Amat et al., 2008). Additionally, many ready-made models 
have been developed, for example, the alignment modules in TxBR 
(Phan et al., 2012), IMOD (Kremer et al., 1996) and our previous 
works (Han et al., 2014). These modules all perform excellently 
in solving the parameter estimation in ET. In this paper, we would 
like to use the alignment module composed in IMOD to illustrate 
the performance of our proposed scheme, because IMOD is a 
successful and widely used suite and is already composed of two 
fiducial marker-based alignment schemes (the half automatic 
detection and tracking scheme appears in version 4.7° and RAPTOR 
of the work by Amat et al. (2008)). However, it should be noted that 
parameter estimation is not the focus of this paper. In practice, read- 
ers can choose alignment modules according to personal preference. 


€ IMOD first automatically detects the fiducial markers in the micrograph with low 
tilt angle based on the user input; IMOD then uses local information to find the 
corresponding markers in the neighboring micrographs; the final tracks should be 
extended and checked by the user. 
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(b) 


Fig. 4. Example images from the test sets. (a) Centriole. (b) Hemocyanin. 


3. Experiments and results 
3.1. Test datasets 


We selected two ET datasets to test the performance of our 
proposed method. 

The first dataset (Centriole dataset) is a tilt series of plastic 
embedded cell section around a centriole region (Fig. 4(a)), which 
was taken on a FEI TF30 microscope (operated at 300 kV) with a 
Gatan Camera. This dataset is obtained from IMOD tutorial’ and 
is an example data to test IMOD fiducial marker based alignment. 
The tilt angles of the projection images range from +65.0° to 
—61.0° at 2° intervals. In total, there are 64 images in the tilt series. 
The size of each tilt image is 1024 x 1024 with a pixel size of 
1.01 nm. The initial orientation of the tilt azimuth with respect to 
the vertical direction of the image is —12.5°. 

The second dataset (Hemocyanin dataset) is a tilt series of vitri- 
fied keyhole limpet hemocyanin solution (Fig. 4(b)). Hemocyanin 
was bought from Sigma-Aldrich (USA) and buffered in PBS solution 
with a protein concentration of 1 mg/ml. The 300 mesh copper EM 
grids with holey carbon film (Quantifoil R2/1) were bought from 
EMS (Electron Microscopy Sciences, USA). Protein solution was 
applied to glow discharged grids and blotted in Vitrobot IV (FEI, 
Netherland) using 4-s blotting time under blotting force 2. The 
blotting chamber conditions were kept at 4 and 100 humidity. 
After blotting, the grids were plunge-frozen in liquid ethane cooled 
by liquid nitrogen. The cryoET data was collected by FEI Titan Krios 
(operated at 300 kV) with a Gatan US4000 camera. The total dose 
used during data collection was kept 8000 e/nm?. There are 95 
images with the tilt angles ranging from —70.0° to +70.0° degree 
at 1° ~ 2° intervals. The size of each tilt image is 2048 x 2048 with 
a pixel size of 0.4nm. The initial orientation of the tilt azimuth 
with respect to the vertical direction of the image is 83.9°. 


3.2. Results 


The first step of our scheme is to refine the value of the marker 
diameter and detect fiducial markers. Fig. 5 shows the average 
template of the fiducial markers of each dataset extracted from 
the slice with 0° tilt angle (the width of the sub-figures is equal 
to the diameter of the markers). We have tested initial input values 
ranging from 0.5 times to 1.5 times the exact value of the marker 
diameter, and the experiment indicated that our method can refine 
the diameter parameter to a suitable value (the value of refined 


7 Downloaded from http://bio3d.colorado.edu/imod/files/tutorialData-1K.tar.gz. 


— 
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Fig. 5. Illustration of the performance when determining the value of marker 
diameter. The sub-figures are average templates of the fiducial markers of the 
Centriole and Hemocyanin datasets. (a) The average marker template of Centriole; 
the diameter value is 10. (b) The average marker template of Hemocyanin; the 
diameter value is 15. 


diameter for Centriole is 9.92 pixels and the value of refined diam- 
eter for Hemocyanin is 15.06 pixels). Fig. 6 illustrates the detected 
fiducial markers from slices with 0° tilt angle and high tilt angle. As 
shown in Fig. 6(a) and (c), our method has detected all of the well- 
distributed markers in the slices with low tilt angles. In Fig. 6(b) 
and (d), the backgrounds of the samples have changed consider- 
ably compared with the backgrounds in Fig. 6(a) and (c); in partic- 
ular, we found that the markers are blurrier in the slice with high 
tilt angle in the Centriole dataset and that the contrast of the back- 
ground has changed drastically in the slice with high tilt angle in 
the Hemocyanin dataset. Because our method considers both the 
“shape” and “contrast” of markers, our method will reject poorly 
shaped the fiducial markers. Thus, our method correctly extracted 
almost all of the well-distributed fiducial markers but rejected the 
poorly-shaped ones, as shown in Fig. 6(b) and (d). Fig. 7 illustrates 
a bar graph of the number of fiducial markers detected by our algo- 
rithm and the number of fiducial markers that we did not detect. 
Notably, our method has a high positive ratio (for most micro- 
graphs, we can detect > 90 percent of the fiducial markers) for 
fiducial marker detection; therefore, the user does not have to pro- 
vide an estimation of the number of fiducial markers as input. Of 
the markers that have not been detected by our method or those 
that were ultimately rejected by our method, these markers can 
be restored during the step of extending the tracks (the method 
described in Section 2.3); this strategy increases the robustness 
of our scheme. 

The second step is to determine the matching of corresponding 
fiducial markers. In our scheme, the matching of micrographs in 
the tilt series obeys the matching strategy developed in 
Section 2.3. However, to demonstrate the robustness of our 
matching method, micrographs separated by large tilt angles 
are used. Fig. 8 illustrates the matching of markers based on the 
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(a) (b) 
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Fig. 6. Examples of detection. All the detected markers are marked by circles, specially, the detected markers that are common to both the 0° and high tilt angle are marked by 
green circles. (a) Detected fiducial markers from the micrograph of Centriole with 0° tilt angle. (b) Detected fiducial markers from the micrograph of Centriole with 64° tilt 
angle. (c) Detected fiducial markers from the micrograph of Hemocyanin with 0° tilt angle. (d) Detected fiducial markers from the micrograph of Hemocyanin with 70° tilt 
angle. 
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Fig. 7. Bar graphs of the number of fiducial markers that were detected and missed by our algorithm. The “bowl” shape of this graph is due to the increased number of fiducial 
markers brought into the field of the image because of foreshortening that caused by the increased tilt angle of the projection. The light blue line shows the viewing area 
enlargement as the tilt angle is increased. (a) Centriole. (b) Hemocyanin. 
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Fig. 8. Example of the matching of corresponding fiducial markers. The Hemocyanin dataset is used here. The moving point set M (labeled by “°”) contains the positions of 
the detected fiducial markers for the micrograph with —55.5° tilt angle. The fixed point set S (labeled by “x”) contains the positions of the detected fiducial markers for the 
micrograph with 9.0° tilt angle. (a) Superimposition of point set M and S. (b) Superimposition of point set T(M) and S, where T(M) means the result dataset of point set M 
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Fig. 9. Histograms of the track lengths for different samples (extracted by our method). At the extreme right hand side of the figures, the large number of tracks represent 
those tracks that cover the entire tilt series. At the left hand side of the figures, there are also a large number of tracks with short track length, which mainly represent the 
points brought in new for the micrographs recorded at higher tilt angle (cannot be tracked throughout the entire tilt series). (a) Centriole. (b) Hemocyanin. 


idea of point set registration. The selected point sets are from the 
Hemocyanin dataset. The moving point set M is comprised of the 
positions of the detected fiducial markers of the micrograph with 
—55.5° tilt angle. The fixed point set S is comprised of the posi- 
tions of the detected fiducial markers of the micrograph with 
9.0° tilt angle. There are totally 134 positions of fiducial markers 
in M, and 84 positions of fiducial markers in S. It is natural that 
more fiducial markers are expected in micrographs with high tilt 
angles. As shown in Fig. 8(a), there is almost no overlap between 
M and S in the superimposition image. Then, M and S are pro- 
cessed by our point set registration algorithm (When considering 
the value of the parameters in Algorithm 4, in our solution, 
dist_thre is set to 5, diameter of wide-base is set to 0.35 of the 
diameter of S and p, is set to 0.999999). In Fig. 8(b), the points 
in M have been transformed to the appropriate positions; an over- 
lap therefore can be observed in the superimposition image of 
T(M) and S. In practice, because of the deformation of the sample 
and the defocus of the micrography, the relationship between 
corresponding markers cannot be described exactly by a simple 
affine transformation; however, for tracking (to obtain the pairs 


of corresponding fiducial markers in different micrographs), the 
deviation is acceptable. In this demonstration, we finally obtained 
74 matching pairs; almost all the pairs of corresponding fiducial 
markers have been found out. 

The length of the tracks is used to demonstrate the overall 
performance of the detection and matching. Fig. 9 illustrates the 
histograms of the track length obtained by our method for the 
two samples (the tilt series were not pre-aligned). In practice, 
our method matches images within one and two view intervals 
and then combines the matching pairs to obtain the tracks. It 
should be noticed that not all of the markers appear throughout 
the entire tilt series; the micrographs with high tilt angles have 
more fiducial markers. However, the fiducial markers in the micro- 
graphs with low tilt angles appear in almost every micrograph of 
the tilt series (regardless of the disappearance caused by the image 
transformation). Therefore, in this work, the ratio between the 
number of tracks with a long length and the number of detected 
fiducial markers in the micrograph with low tilt (zero tilt) is 
used to indicate the effectiveness of our detection and tracking. 
As shown in Fig. 9, there are 38 tracks covering more than 
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Fig. 10. Histograms of the track lengths for different samples (extracted by RAPTOR). (a) Centriole. (b) Hemocyanin. 


54 micrographs (approximately 85% of the entire tilt series) in the 
Centriole dataset, and 69 tracks covering more than 81 micrographs 
(approximately 85% of the entire tilt series) in the Hemocyanin 
dataset. Furthermore, we obtained 27 tracks (64 micrographs in 
length) and 60 tracks (95 micrographs in length) that cover the 
entire tilt series of the Centriole and Hemocyanin datasets, 
respectively. As a comparison, the numbers of detected fiducial 
markers in the micrographs of the Centriole and Hemocyanin 
datasets with zero tilt angle are 51 and 84, respectively. 
Consequently, the tracks with the entire length have covered 
74.5 percent (38/51) and 82.1 percent (69/84) of the fiducial 
markers in the zero tilt micrographs of the Centriole and 
Hemocyanin datasets, respectively. If excluding the disappearance 
of fiducial markers nearing the border of the micrographs (caused 
by image transformation), the value of coverage can be higher. The 
high value of coverage indicates the effectiveness of our method. 

For comparison, we also analyzed the samples using the half 
automatic marker detection and tracking scheme provided by 
IMOD 4.7 and the scheme of RAPTOR (Amat et al., 2008). We 
compared the number of tracks that covered the entire tilt series 
and were generated by our method, IMOD and RAPTOR. In the 
half automatic marker detection and tracking scheme of IMOD, 
31 and 55 tracks of fiducial markers can be tracked for the 
Centriole and Hemocyanin datasets, respectively, all of which 
cover the entire tilt series. Thus, compared with the result of 
IMOD, our scheme does not perform as well for the Centriole 
dataset (27 vs. 31) but performs better for the Hemocyanin 
dataset (60 vs. 55). However, to obtain these tracks, our operator 
has interacted with IMOD a dozen times to relocate and retrack 
the markers (half automatic); in our scheme, no manual interac- 
tion is required. In the scheme of RAPTOR, 11 and 46 tracks of 
fiducial markers can be tracked for the Centriole and Hemocyanin 
datasets, respectively, all of which cover the entire tilt series. 
Further statistics of the length of the tracks detected by RAPTOR 
are shown in Fig. 10 (a total of 39 and 77 tracks for the Centriole 
and Hemocyanin datasets, respectively). For the Centriole dataset, 
the number of tracks with the entire length detected by our 
method is 27, while the number of tracks with the entire length 
detected by RAPTOR is 11. For the Hemocyanin dataset, the 
number of tracks with the entire length detected by our method 
is 60, while the number of tracks with the entire length detected 
by RAPTOR is 46. Therefore, the quality of tracks detected by our 
method is much higher than the ones detected by RAPTOR 
(as shown in Figs. 9 and 10, the total number of tracks detected 
by our method is also larger than the one of RAPTOR). 


The final step is to calculate projection parameters. Here, we 
use the alignment module provided by IMOD (Kremer et al., 
1996). In the bundle adjustment, our scheme only uses the tracks 
that cover more than 85 percent of the micrographs (thus, all the 
used markers have corresponding markers in the zero tilt 
micrograph) to calibrate the projection parameters. Fig. 11 
illustrates the physical trajectories of all selected tracks in image 
space (x-y coordinates in pixel number). No pre-alignment was 
performed on the datasets. Fig. 11(a) and (c) indicate that there 
are substantial initial rotations in both datasets and that the 
images in the tilt series vibrate severely. Despite the diverse 
appearance of these two datasets, our scheme can accurately 
detect and track the fiducial markers. The results illustrated in 
Fig. 11(b) and (d) demonstrate that the alignment based on our 
final tracks is successful. Reprojection residuals can be used to 
evaluate the quality of the extracted fiducial markers. For 
comparison, we also analyzed the samples using RAPTOR 
(Amat et al., 2008) and half automatic marker selection and 
tracking scheme in IMOD. Detailed information about the 
alignment results are summarized in Table 1. For the Centriole 
dataset, our scheme used 38 tracks as the input of the bundle 
adjustment; the scheme of IMOD used 31 as the input of bundle 
adjustment and the scheme of RAPTOR used 39 tracks (the 
analysis of tracks has been illustrated in previous paragraph). 
As shown in Table 1, for the Centriole dataset, the initial residual 
of our method is close to the results of IMOD and better than the 
results of RAPTOR. The residual of Ali. of our method is even 
smaller than the result of IMOD. The results of our method and 
IMOD are much better than the results of RAPTOR. For the 
Hemocyanin dataset, our scheme used 69 tracks as the input of 
the bundle adjustment; the scheme of IMOD used 55 as the input 
of bundle adjustment and the scheme of RAPTOR used 77 tracks. 
As shown in Table 1, for the Hemocyanin dataset, the overall 
results of our method is almost identical to the results of IMOD 
and much better than the results of RAPTOR. Fig. 12 illustrated 
the reconstruction result of each sample. Judging from Fig. 12, 
the alignment and reconstruction results carried out by our 
scheme are as good as the results obtained by manual detection 
and tracking of an expert in details. However, the alignment 
result of our scheme is obtained under the situations fully 
automatic and totally without any manual interference. There- 
fore, we can draw the conclusion that our fully automatic scheme 
is much easier to use and can achieve results comparable with 
IMOD, much better than RAPTOR. 
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Fig. 11. (a) Overlay of raw fiducial marker tracks in image space (x-y coordinates in pixel number) extracted from the Centriole dataset. (b) Overlay of the aligned fiducial 
marker tracks of the Centriole dataset in image space. (c) Overlay of raw fiducial marker tracks in image space extracted from the Hemocyanin dataset. (d) Overlay of the 


aligned fiducial marker tracks of the Hemocyanin dataset in image space. 


Table 1 
Summary of residuals carried out by different methods. 
Method* Our scheme RAPTOR IMOD 
Ali. Ali.* Ali. Ali.4 Ali. Ali.* 
Residual” avg. std. avg. std. avg. std. avg. std. avg. std. avg. std. 
Dataset: Centriole 0.57 0.42 0.46 0.26 0.75 0.51 0.68 0.38 0.58 0.38 0.51 0.30 
Dataset: Hemocyanin 0.71 0.51 0.63 0.47 0.86 0.56 0.77 0.47 0.70 0.48 0.63 0.44 


è The alignment module provided by IMOD can fix fiducial markers with large residuals through manual adjustment. Here, “Ali. 
of bundle adjustment, i.e., the initial alignment (there are no manual adjustment for large residual markers in initial alignment); “Ali. 


1” refers to the alignment result after a turn 


4» refers to the alignment result after 4 


turn of adjustment (a turn of alignment means a turn of manual adjustment for large residual markers and a turn of bundle adjustment). 
> The abbreviation “avg.” refers to the mean value of residuals; “std.” refers to the standard deviation of residuals. 


4. Discussion and conclusion 


In this paper, we present a novel, fully automatic scheme for 
alignment in ET. Our scheme focuses on the detection and tracking 
of fiducial markers. In our scheme, the refining of the marker 
diameter parameter and the idea that further converting the 
detection problem to sampling and classification problem provide 
an accurate detection of fiducial markers with a high recognition 
ratio. The procedure of refining fiducial marker localization pro- 
vides highly precise positions of fiducial markers. The matching 
of corresponding markers based on the idea of incomplete point 
set registration under affine transformation guarantees a global 
matching and tracking of fiducial markers. Furthermore, the reduc- 
tion of the matching problem to a point set registration problem 


could provide a novel solution for other similar problems in the 
field of structural biology. The experimental results indicate that 
our method can provide excellent tracking quality, with almost 
identical results to the half automatic scheme of IMOD and better 
results than RAPTOR. Our scheme is fully automatic, depends on 
fewer parameters (only requires the marker diameter to be 
inputted by the user) and does not require any manual interaction. 
Therefore, the proposed scheme is much easier to use than the half 
automatic scheme of IMOD and provides the possibility of auto- 
matic batch processing of electron tomographic reconstruction. 
Along with the new techniques introduced for the processing of 
electron tomography data, the most important contribution of our 
work to the scope of ET is its “fully automatic” nature, which can 
completely free the researchers from tedious data processing. 
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Fig. 12. (a) The central x-y section of the reconstruction of the Centriole dataset aligned by our scheme. (b) A similar x-y section of the Centriole dataset aligned by the half 
automatic scheme of IMOD. (c) The central x-y section of the reconstruction of the Hemocyanin dataset aligned by our scheme. (d) A similar x-y section of the Hemocyanin 
dataset aligned by the half automatic scheme of IMOD. The sharpness of details is similar in both reconstructions. 


Especially, concerning the data acquisition in subtomogram aver- 
aging (which requires handling dozens of tilt series with the iden- 
tical protein or ultrastructure), our scheme can save researchers a 
substantial amount of time. On the other hand, for tilt series with a 
complex projection model (for example, the non-linear bundle 
adjustment of Lawrence et al. (2006) requires a large number of 
fiducial markers), our scheme can also provide enough coverage 
of fiducial marker detection and good fiducial marker tracks in a 
fully automatic way. 

Although the effectiveness of our method has been shown, sev- 
eral points must be addressed. 

First, the appropriate affine transformation between two views 
of tilt series is assumed. Although the appropriate affine transfor- 
mation relationship of two respective 2D views from the identical 
planar surface has been widely used in computer vision 
(Huttenlocher and Ullman, 1990; Koenderink and van Doorn, 
1991), it may not be able to describe relatively large deforma- 
tions, mass losses or other distortion effects of samples. Thus, 
the assumption may limit the application of our method. In the 
situation with large distortions, the results of our method can 
be used as the initial input of the iterative closest point (ICP) 
with more parameterized models (Besl and McKay, 1992; Jian 
and Vemuri, 2011). However, in almost all experiments, the 


micrographs with large distortions may directly lead to a failure 
of 3D tomography. Therefore, in most situations, our method does 
not encounter this problem. 

Second, although some techniques based on local informations 
can be introduced to the model of point set registration and 
accelerate the computational speed, these techniques are not 
recommended because they may introduce limitations to our 
application. For example, technique of feature descriptors (e.g. 
Scale-Invariant Feature Transform (SIFT)) can reduce the execu- 
tion time of our method in samples with sufficient biological 
structures, the hybrid scheme will not be able to applied in the 
condition that the biological structures are insufficient. If we 
introduce local informations (e.g. an estimated corresponding 
search area from the dataset after coarse alignment), the roughly 
inaccurate matching pairs obtained by local search could provide 
as the seeds for estimating the transformation matrices in the 
point set registration, consequently, we must bear the risk of 
the failure of local search, which will hamper the robustness of 
our scheme. 

Finally, our algorithms proposed here may be more feasible in 
the alignment schemes of dual-axis or multi-axis electron 
tomography. Our method can smoothly solve the fiducial marker 
tracking problem across different rotation axes. In other words, 
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our method can naturally determine the corresponding relation- 
ship of fiducial markers across different micrographs, in a globally 
consistent manner. This capability can make it easy to perform 
bundle adjustment in a unique reference system, or make it 
easier to simultaneously align a dual-axis tilt series based on ref- 
erence points (Cantele et al., 2010). 


5. Implementation 


Software that is compatible with IMOD and based on our 
method, markerauto (fully automatic marker-based alignment), is 
being developed. Interested users can download the beta version 
in our site http://ear.ict.ac.cn/. 
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Appendix. Algorithms’ flowcharts 
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Fig. A.1. Flowchart of the marker detection procedure with the assumption that the 
value of marker diameter is already known. 
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Fig. A.2. Flowchart of a loop of the marker diameter refinement procedure with the 
assumption that the detected fiducial markers are real fiducial markers and 
accurately localized. 
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Fig. A.3. Flowchart of the marker position refinement procedure. 
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